clear all;
clc;
%%
% Hochrechnung auf eine Sekunde Bestrahlung (Soll: 10^16  Photonen)
simfaktor = 1e9;
nDiv = 400;
nDivX = nDiv;
nDivY = nDiv;
nDivZ = nDiv;
view = 'front';

% Neu!!!!
energScale = [1, 1, 1];
% energScale = [20/3, 3, 1];

% Position Mittelpunkt Tumor
tumorX = ceil(0.5*nDivX);
tumorY = ceil(250/400*nDivY);
tumorZ = nDivZ-ceil(190/400*nDivZ);

tumorRadius = 20; % mm

% Laenge eines Voxels in mm
lenX = 400/nDivX;
lenY = 400/nDivY;
lenZ = 400/nDivZ;

switch(view)
    case 'front'
        schnitt1 = tumorY;
        schnitt2 = tumorX;
        nDivT = nDivZ;
        tumorT = tumorZ;
        lenT = lenZ;
    case 'top'
        schnitt1 = tumorX;
        schnitt2 = tumorZ;
        nDivT = nDivY;
        tumorT = tumorY;
        lenT = lenY;
    case 'side'
        schnitt1 = tumorY;
        schnitt2 = tumorZ;
        nDivT = nDivX;
        tumorT = tumorX;
        lenT = lenX;
end

if(exist('filename', 'var') == 0)
    [filename, pathname] = uigetfile('.log', 'MultiSelect', 'on');
    
    if iscell(filename)
        kEnd = numel(filename);
    else
        kEnd = 1;
    end
    
    dat = zeros(nDiv, kEnd); %TODO
    
    for k=1:kEnd
        tempDat = zeros(nDiv, 1);
        
        if iscell(filename)
            file = fopen(fullfile(pathname, filename{k}), 'r');
        else
            file = fopen(fullfile(pathname, filename), 'r');
        end
        
        % Datedatei mit Kopf-Daten
        % Name (xCoord*nDivY*nDivZ+yCoord*nDivZ+zCoord) eDeposit doseDeposit
        input = fscanf(file,'%u\t%e\t%e\n',[3 inf]);
        
        % Bitte daran denken das File zu schliessen, sonst gibt es Probleme
        fclose(file);
        
        % Anzahl Eintraege/"Zeilen" in aktuell eingelesener Datei
        entries = numel(input(1,:));
        
        for n=1:entries
            nx = floor(input(1,n) / (nDivY*nDivZ));
            ny = floor((input(1,n) - nx*nDivY*nDivZ) / nDivZ);
            nz = input(1,n) - nx*nDivY*nDivZ - ny*nDivZ;
            
            switch(view)
                case 'front'
                    indexSchnitt1 = ny+1;
                    indexSchnitt2 = nx+1;
                    indexT = nz+1;
                case 'top'
                    indexSchnitt1 = nx+1;
                    indexSchnitt2 = nz+1;
                    indexT = ny+1;
                case 'side'
                    indexSchnitt1 = ny+1;
                    indexSchnitt2 = nz+1;
                    indexT = nx+1;
            end
            
            if ((indexSchnitt1 == schnitt1) && (indexSchnitt2 == schnitt2))
                tempDat(indexT) = input(3,n)*simfaktor;
%                 dat(indexT,k) = input(3, n) * simfaktor;
            end
        end
        
        dat(:, k) = tempDat * energScale(k);
        
        clear input
    end
    
end

space = [0:lenZ/10:(nDivZ-1)*lenZ/10];

maxVal = max(max(dat));

%% Tumor im Plot mar
plotTumorLeft = (tumorT - tumorRadius/lenT)/10;
plotTumorRight = (tumorT + tumorRadius/lenT)/10;
plotTumorTop = 0;
plotTumorDown = maxVal;
plotHeadLeft = 5;
plotHeadRight = 35;
plotHeadTop = plotTumorTop;
plotHeadDown = plotTumorDown;


%% Plot
figure;
axes('FontSize', 14);
plot(space,dat)

hold on
plot([plotTumorLeft, plotTumorLeft], [plotTumorDown, plotTumorTop],'k','Marker','.','LineWidth',3)
plot([plotTumorRight, plotTumorRight], [plotTumorDown, plotTumorTop],'k','Marker','.','LineWidth',3)
plot([plotHeadLeft, plotHeadLeft], [plotHeadDown, plotHeadTop],'--k','Marker','.','LineWidth',3)
plot([plotHeadRight, plotHeadRight], [plotHeadDown, plotHeadTop],'--k','Marker','.','LineWidth',3)
hold off

axis([0 40 0 0.02])

%annotation('textbox',[tumorZ/10 maxVal/2 10 maxVal/4],'String','Tumor')
annotation('textbox',[0.47 0.535 0.13 0.065],'String','Tumor', 'HorizontalAlignment', 'center', 'LineStyle', 'non', 'FontSize', 14)
if strcmp(view, 'front') || strcmp(view, 'side')
    annotation('textbox',[0.16 0.535 0.23 0.065],'String','Head Front', 'HorizontalAlignment', 'center', 'LineStyle', 'non', 'FontSize', 14)
    annotation('textbox',[0.39 0.535 0.7 0.065],'String','Head Back', 'HorizontalAlignment', 'center', 'LineStyle', 'non', 'FontSize', 14)
end

legend(filename,'Location','NorthEast');
%axis([0 40 0 100]);
xlabel('Entfernung von der Strahlenquelle in cm');
ylabel('Energieabgabe in Gy');
title('Vergleich: Energieabgabe der Stahlung in Abhaengigkeit der Tiefe');

if kEnd == 1
    maxVorTumor = max(dat(1:ceil(tumorT - tumorRadius/lenT)))
    maxInTumor = max(dat(ceil(tumorT - tumorRadius/lenT):ceil(tumorT + tumorRadius/lenT)))
    minInTumor = min(dat(ceil(tumorT - tumorRadius/lenT):ceil(tumorT + tumorRadius/lenT)))
    maxHinterTumor = max(dat(ceil(tumorT + tumorRadius/lenT):nDivT))
end

dx = lenZ/10;
KEin = 5/dx; KAus = 35/dx; TEin = plotTumorLeft/dx; TAus = plotTumorRight/dx;

EnergieAbgabe = zeros(4,3);
EnergieAbgabe(1,:) = sum(dat(KEin:TEin,:));
EnergieAbgabe(2,:) = sum(dat(TEin:TAus,:));
EnergieAbgabe(3,:) = sum(dat(TAus:KAus,:));
EnergieAbgabe(4,:) = sum(dat(KEin:KAus,:));
AnteilDosisInTumorZuDosisKopf = EnergieAbgabe(2,:) ./ EnergieAbgabe(4,:);
